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Abstract 

We have performed Monte Carlo simulations of the Ising model coupled to three- 
dimensional quantum gravity based on a summation over dynamical triangulations. 
These were done both in the microcanonical ensemble, with the number of points 
in the triangulation and the number of Ising spins fixed, and in the grand canonical 
ensemble. We have investigated the two possible cases of the spins living on the vertices 
of the triangulation ( "direct" case) and the spins living in the middle of the tetrahedra 
("dual" case). We observed phase transitions which are probably second order, and 
found that the dual implementation more effectively couples the spins to the quantum 
gravity. 







1 Introduction 



Following recent Monte Carlo simulations of spin models coupled to two-dimensional quan- 
tum gravity [|I|] and of pure three-dimensional quantum gravity we have performed some 
small exploratory simulations of the Ising model coupled to three-dimensional quantum grav- 
ity. This is the simplest model of matter coupled to quantum gravity in three dimensions, 
and provides a step on the way to simulating the physical case of matter coupled to four- 
dimensional quantum gravity. 

The continuum Einstein action for 3-dimensional gravity is 

s^-l^'iVs{x-^), (1) 

where A is the cosmological constant and G is Newton's constant. When a dynamical trian- 
gulation (with Nq vertices, A^i links, N2 triangles and N^, tetrahedra) is used to regularize 
the functional integration over metrics the two terms in (|I]) discretize as follows: 

Jd'^^^Ns (2) 

/ d3e^i?^5:(c-n(/)), (3) 

I 

where n{l) is the number of tetrahedra which share the link /, and c = 27r/a, with a = 
arccos{^), is a number which ensures that R is zero for flat space. Since each tetrahedron 
is comprised of six links, X); '^(O — ^-^3 ^"^^ (0) can be written as cNi — Hence the 

discrete version of (|5) is 

^ = AiVg - ^ (ciVi - 6Ns) . (4) 

Now since we have the following relations for a three-dimensional triangulation (or "three- 
dimensional simplicial manifold") 

N3-N2 + Ni-No = (5) 

N2 = 2N3 (6) 

(first is Euler's relation; second is because each triangle is common to two tetrahedra), we 
can change to the more convenient variables Nq, and write 

S = AaiVg - AoiVo, (7) 

with A3 = (A + 1^) and Aq = Therefore the grand canonical partition function for 

quantum gravity in three dimensions based on a summation over dynamical triangulations 
T with tetrahedra and Nq vertices is 



^(Aa, Ao) = J dNsdNoZN^^N^e' 
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where 




E PiT) 



(9) 



T(N3,No} 



is the micro-canonical partition function in which N^jNq are fixed, and p{T) is a suitable 
weight for each triangulation T. One usually appeals to universality and sets p(T) = 1. 

In order to perform a Monte Carlo simulation of either of these partition functions for 
pure three-dimensional quantum gravity one must define a set of "moves" which ergodically 
sum over all possible triangulations. In two dimensions this is straight-forward since one 
move suffices for each case: for the grand canonical partition function one uses the so-called 
"split-join" move, and for the micro-canonical partition function the "fiip". Unfortunately 
things are not as simple in three dimensions - there is no known set of moves which is ergodic 
in the micro-canonical ensemble; however two moves exist for the grand canonical case: the 
so-called "Alexander move" which is a generalization of the two-dimensional split-join move 
changing one tetrahedron into four tetrahedra (or vice versa), and a generalization of the 
'fiip" taking two tetrahedra into three tetrahedra (or vice versa). For more details and 
illustrations of these moves see 0. 

If we now couple the Ising model to this discretization of three-dimensional quantum 
gravity the action (0) becomes 



where Sij is the connectivity matrix of the triangulation, (Xj are the Ising spins taking values 
-|-1 or —1, and (3 is the inverse temperature. There are two possible ways in which to add the 
Ising spins to the dynamical triangulation: they can either live on the vertices in which case 
there will be A^o spins or they can live in the tetrahedra giving spins. We shall refer to 
the former as the "direct" case and the latter as the "dual" . For the direct implementation 
there are Ni = Nq + N3 links between the spins so Sij will contain this number of unit 
entries; for the dual case there are N2 = 2Ns interactions, and hence non-zeros in Sij. Note 
that we have arranged the normalization such that at zero temperature, i.e. f3 = 00, all the 
Ising spins line up and their contribution to Sj vanishes. At /3 = we recover pure quantum 
gravity. 

As for the pure three-dimensional quantum gravity case, there are two possible partition 
functions: grand canonical and micro-canonical. For the grand canonical ensemble, one 
must create and destroy Ising spins when either vertices in the direct case or tetrahedra 
in the dual case are created and destroyed. Destruction of Ising spins is trivial - they are 
simply removed - however creation is rather more subtle - the values of the new spins 
should be picked randomly to avoid any bias. Our simulation of the grand canonical Ising 
plus three-dimensional quantum gravity uses the same methods as were used in all the pure 
three-dimensional quantum gravity simulations, see |^ for details. 

As a step on the way to the grand canonical simulation, we have chosen to simulate 
the micro-canonical ensemble. To do this N^, Nq and the number of Ising spins must be 
fixed. Therefore we use only the fiip move to update the triangulation, and we use it in 
pair-wise fashion - firstly changing two tetrahedra to three tetrahedra then immediately 




(10) 
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after (somewhere else in the triangulation) changing three tetrahedra to two tetrahedra. As 
explained above this is not known to be ergodic for the pure three-dimensional quantum 
gravity case, so we cannot assume it to be when we couple in Ising spins. However it is 
an easier simulation to do and may provide an approximation to the more realistic grand 
canonical case (this is somewhat like simulating quenched QCD as an approximation to full 
QCD in lattice gauge theory). 

For both the grand canonical and micro-canonical simulations, in addition to updating the 
triangulation, we must update the Ising spins. To do this we make use of cluster algorithms, 
introduced in 0, which are very effective in the vicinity of phase transitions, where spins 
tend to form large clusters. However, away from the transition, the standard Metropolis 
method is more practical. Therefore, we used both algorithms to update the Ising spins, 
doing one pass of each per update step. We actually used Wolff's version of the cluster 
algorithm as it is faster than Swendsen- Wang's [^. 



2 Micro- canonical 

Physically the vacuum in four-dimensional gravity is almost flat, therefore in our micro- 
canonical simulations of three-dimensional quantum gravity we pick N3, No so that the av- 
erage curvature of the triangulation 

27r(iV3 + iVo)-6aiV3 

is as small as possible; which in turn implies that ^ ~ .1755. This leads to the following 
values of A^o, N^. 69,393 ; 154,878 and 325,1852. If is the number of tetrahedra meeting 
at the vertex i then, since each tetrahedron has four vertices, J2i^3 = 4:^3^ and therefore 
< >= 4^ ^ 22.8. Thus each vertex is shared by roughly 23 tetrahedra. This means 
that for the direct implementation of the Ising spins on the triangulation, each spin will have 
on average 23 neighbors, whereas for the dual implementation each spin will have only four 
neighbors (the number of faces of a tetrahedron). Therefore we can expect a much stronger 
correlation of spins for a given (3 for the direct case compared to the dual case, which causes 
the significantly stronger phase transition to occur at a lower critical value of (3. We shall 
first present our results for the direct case and then for the dual case. In both cases we 
measured the standard thermodynamic energy E and magnetization M for the Ising model, 
the acceptance rate of the Metropolis algorithm qm, the acceptance rate for the flip moves 
ttf, and the fractal dimension of the spin clusters constructed by the Wolff algorithm. 

For the direct case the number of Ising spins is A^^q = 69, 154 and 325. We ran simulations 
at the following 18 values of f3: 0.01, 0.02, 0.03, 0.0325, 0.035, 0.0375, 0.04, 0.0425, 0.045, 0.0475, 
0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.2, 0.3. The results for E, M, om, aj for the A'o = 325 simula- 
tion are listed in Table 1. 
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Table 1: Measured values of energy E, magnetization M, Metropolis acceptance au and flip 
acceptance a/ from Nq — 325 direct simulation. 

By numerically differentiating the energy and magnetization we obtain the specific heat 
and susceptibility shown for all A^o in Figs. 1 and 2 respectively. Using a cubic spline 
approximation we estimate the positions Cmax and Xmax of the maximum of the peaks to 
obtain the critical inverse temperatures /3c for each A^o hsted in Table 2. 
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Table 2: Position of peaks in specific heat and susceptibility for each value of A^o from direct 
simulation. 

We could attempt to fit to the specific heat and susceptibihty peaks in order to extract 
exponents and/or we could try finite-size scaling analysis but as these are exploratory sim- 
ulations of small systems the results would not be very impressive or reliable. Therefore we 
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leave this for future calculations. Here we can say only that there is a phase transition, at 
approximately (3c = 0.048(2) for Nq = 325, which is probably second order since the peaks 
in C and x increase with system size (first order is ruled out as there are no discontinuites 
in E and M). Unfortunately without the critical exponents it is not possible to extrapolate 
the (3c values to A^o = oo. 

Turning to the acceptance rates we see that the Metropolis acceptance falls off as (3 
increases as expected - the spins become frozen as the temperature falls. However the flip 
acceptance rate does not change very much from its initial value of 1 at /3 = (where the 
spins are random so do not systematically effect the triangulation) . In other words the Ising 
spins do not seem to interact very strongly with three-dimensional quantum gravity when 
they are implemented directly - i.e. live on the vertices of the triangulation. Numerically 
this is obviously due to the fact that during the flip move only one out of on average 23 Ising 
spins at each vertex is reconnected so the energy change is relatively small. 

Now we turn to the results for the dual case, where the number of Ising spins is N-^ = 
393,878 and 1852. Simulations were done at 12 values of (3: 0.1,0.2,0.3,0.325,0.35,0.375, 
0.4, 0.425, 0.45, 0.475, 0.5. 0.6: the results from the .¥3 = 878 siniTilation are in Table 3. 
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Table 3: Measured values of energy E, magnetization M, Metropolis acceptance qm and flip 
acceptance a/ from = 878 dual simulation. 

We show the specific heat and susceptibility, obtained from numerical differentiation, for 
A^3 = 393,878 (we do not have enough data for = 1852), in Figs. 3 and 4 respectively. 
And we estimate the position of the maximum of the peaks - Table 4. 
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Table 4: Position of peaks in specific heat and susceptibility for two values of N3 from dual 
simulation. 

Again it is not worth attempting to fit to this data to extract critical exponents. For now we 
can say that we have seen a second order phase transition at approximately jSc = 0.418(2) 
for N3 = 878. 

Thus for the dual case, the Ising spins significantly affect the triangulation - the accep- 
tance rate for the flip move drops from near 1 at low (3 to almost as /? increases above 
f3c- Thus matter in the form of Ising spins interacts strongly with quantum gravity when 
it is coupled to it through the dual of the triangulation - i.e. when the spins live in the 
tetrahedra. If we compare the direct and dual cases we see that for the former the peaks 
in C and x ci-re much larger as expected. Interestingly the critical value of /3 for the dual 
case in which each spin is connected to four neighbors {jSc = 0.418(2)) is closer to the value 
(0.221652(3) [^) for the standard 3D Ising model which has coordination number six, than 
is the value (0.048(2)) from the direct case with its average coordination number of 23. 

Finally, we look at the fractal dimension of the clusters constructed by the Wolff algo- 
rithm. This is interesting because, as discussed by Meo, Heermann and Binder ||^, these 
clusters are actually the real, physical clusters in the system. For a standard Ising model 
simulated in D dimensions, the fractal dimension d is given hy d = D — f3/v which is approx- 
imately 2.5 for D = 3 0. In the simulations we stored the sizes and diameters Q of the Ising 
clusters. The log — log plots of the average cluster size versus it's diameter along with the 
best fits for d, for both the direct and dual simulations, are shown in Fig. 5. For the direct 
case we get a value oi d = 2.7(2) close to the standard 3D Ising model value, for the dual 
case it is a little smaller (2.0(1)). These values should be considered as preliminary given 
the small triangulations used - the maximum sized clusters which fit on the direct and dual 
triangulations are only about 50 and 100 sites respectively. 

3 Grand canonical 

Armed with the knowledge that the dual implementation of the Ising spins on the trian- 
gulation more effectively couples them to three-dimensional quantum gravity in the micro- 
canonical ensemble, we performed a modest simulation of the dual case in the grand canon- 
ical ensemble. To do this one first picks likely values of the parameters A3,Ao plus widths 
AA3, AAo then lets the Monte Carlo simulation vary the former within the latter and find 
the actual values of A3, Aq. We started with A3 = l,Ao = 0,AA3 = AAq = 0.25 and ran 
the first simulation with N^^ = 393 at p = 0.1,0.2,0.3,0.35,0.4,0.45,0.5. This resulted in 
the predictions for A3, Aq listed in columns 2 and 3 respectively in Table 5. Then we used 
these predictions as start values for the second simulation with N3 = 878 resulting in the 
further predictions listed in columns 4 and 5. Also listed in Table 5, in the "/? = 0" row, 
are the values of A3, Aq from simulations of pure three-dimensional quantum gravity on the 
same-sized triangulations 0. In these pure quantum gravity simulations a phase transition 

^By diameter we mean the distance from the random starting point used in construction of the cluster to 
the most remote point on the cluster border; on average this should be proportional to the gyration radius. 



6 



was found at Aq ~ 0.61; for Aq < 0.61 we are in a "strong gravity" phase (Newton's constant 
G ~ ^ is large) where the fractal or Hausdorff dimension of the triangulation becomes 
large. (Note dn should not be confused with the fractal dimension of the Ising spin clus- 
ters.) In Fig. 6 we see a roughly linear relationship between the A3,Ao values and (3 with 
no dramatic change at the Ising model phase transition around I3c = 0.4. The Ising model 
reduces the values of both A3, Aq and therefore has the effect of reducing the cosmological 
constant and making already strong gravity stronger. Perhaps this makes sense - matter 
does couple positively to gravity after all, and the cosmological constant is essentially zero 
in the real universe. However, it is also apparent that there are large finite-size corrections 
which tend to increase A3 and decrease Aq. Therefore when large simulations are done it may 
turn out that Aq remains at its pure quantum gravity value despite the addition of matter 
in the form of Ising spins, whereas A3 is significantly reduced, perhaps to zero. 
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Table 5: Predicted values of A3, Aq from grand canonical dual simulations with = 393 
and 878. 

In addition to E, M, aM,ctf we also measured the average curvature of the triangulation. 
In the strong gravity regime for pure quantum gravity this is found to be negative [@ and 
indeed we find this here. The values of E, M and qm for the grand canonical simulation are 
consistent with those for the micro-canonical; therefore, despite doubts about its ergodicity, 
the micro-canonical simulation does provide a good approximation to the grand canonical. 
The values of a/ are substantially reduced because now this quantity includes the acceptance 
rate for the Alexander moves as well as the flip moves, and the former are greatly suppressed 
by the relatively large change in energy brought about by the creation and destruction of 
(in this case three) Ising spins. 

The fractal dimension of the clusters constructed by the Wolff algorithm, remains close 
to that found for the micro-canonical simulation, namely a little less that the value for the 
standard 3D Ising model. This is probably due to the fact that the underlying triangulation 
has a large Hausdorff dimension, therefore clusters of spins built upon it grow more slowly 
than those built with an underlying 3D lattice. 
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Table 6: Measured values of energy E, magnetization M, Metropolis acceptance om, flip 
acceptance a/ and average curvature R from A^3 = 878 grand canonical dual simulation. 

We have not yet investigated the more physically relevant case of "weak gravity" where 
Newton's constant is small and Aq large (at least > 0.61). This would involve fixing Aq and 
hence A^o, then performing a canonical simulation in which A3 and A'3 are allowed to vary. 
Here we have instead let the system pick its preferred value of Aq. 

4 Conclusions 

Wc have investigated the properties of the Ising model coupled to three-dimensional quan- 
tum gravity. Firstly, from micro-canonical simulations, we have determined that the dual 
implementation of Ising spins in the tetrahedra of the triangulation is more effective than the 
direct implementation on vertices, and that both cases exhibit a phase transition which is 
most likely second order. Then, for the dual implementation, we performed a grand canoni- 
cal simulation which also displays a second order phase transition, and discovered that the 
effect of the Ising spins is to definitely reduce the cosmological constant and perhaps make 
already strong quantum gravity stronger. 

The next step in these simulations is to exhaustively explore the parameter space by 
choosing fixed values for Aq both in the weak and strong pure quantum gravity phases, and 
- of course - to simulate larger systems so that finite-size scaling analysis can be done to 
obtain critical exponents and extrapolate to the continuum limit. 
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Figure Captions 



Fig. 1. Specific heat for direct micro-canonical simulation. 

Fig. 2. Susceptibility for direct micro-canonical simulation. 

Fig. 3. Specific heat for dual micro-canonical and grand canonical simulations. 

Fig. 4. Susceptibility for dual micro-canonical and grand canonical simulations. 

Fig. 5. Cluster size vs. diameter, at closest j3 values to f3c, for direct No = 325 and dual 
A^3 = 878 micro-canonical simulations, along with fits to d shown as dotted and full 
lines respectively. 

Fig. 6. Predictions for Aq, A3 from grand canonical simulations. 
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